home *** CD-ROM | disk | FTP | other *** search
/ Aminet 35 / Aminet 35 (2000)(Schatztruhe)[!][Feb 2000].iso / Aminet / gfx / misc / gnuplot-src.lha / gnuplot-3.7.1src / gnuplot-3.7.1.lha / gnuplot-3.7.1 / fit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-10-01  |  45.7 KB  |  1,650 lines
  1. #ifndef linti < num_params; Sid = "$Id: fit.c,v 1.23.2.3 1999/10/01 18:22:14 lhecking Exp $";
  2. #endif
  3.  
  4. /*
  5.  *    Nonlinear least squares fit according to the
  6.  *      Marquardt-Levenberg-algorithm
  7.  *
  8.  *      added as Patch to Gnuplot (v3.2 and higher)
  9.  *      by Carsten Grammes
  10.  *      Experimental Physics, University of Saarbruecken, Germany
  11.  *
  12.  *      Internet address: cagr@rz.uni-sb.de
  13.  *
  14.  *      Copyright of this module:  1993, 1998  Carsten Grammes
  15.  *
  16.  *      Permission to use, copy, and distribute this software and its
  17.  *      documentation for any purpose with or without fee is hereby granted,
  18.  *      provided that the above copyright notice appear in all copies and
  19.  *      that both that copyright notice and this permission notice appear
  20.  *      in supporting documentation.
  21.  *
  22.  *      This software is provided "as is" without express or implied warranty.
  23.  *
  24.  *      930726:     Recoding of the Unix-like raw console I/O routines by:
  25.  *                  Michele Marziani (marziani@ferrara.infn.it)
  26.  * drd: start unitialised variables at 1 rather than NEARLY_ZERO
  27.  *  (fit is more likely to converge if started from 1 than 1e-30 ?)
  28.  *
  29.  * HBB (Broeker@physik.rwth-aachen.de) : fit didn't calculate the errors
  30.  * in the 'physically correct' (:-) way, if a third data column containing
  31.  * the errors (or 'uncertainties') of the input data was given. I think
  32.  * I've fixed that, but I'm not sure I really understood the M-L-algo well
  33.  * enough to get it right. I deduced my change from the final steps of the
  34.  * equivalent algorithm for the linear case, which is much easier to
  35.  * understand. (I also made some minor, mostly cosmetic changes)
  36.  *
  37.  * HBB (again): added error checking for negative covar[i][i] values and
  38.  * for too many parameters being specified.
  39.  *
  40.  * drd: allow 3d fitting. Data value is now called fit_z internally,
  41.  * ie a 2d fit is z vs x, and a 3d fit is z vs x and y.
  42.  *
  43.  * Lars Hecking : review update command, for VMS in particular, where
  44.  * it is not necessary to rename the old file.
  45.  *
  46.  * HBB, 971023: lifted fixed limit on number of datapoints, and number
  47.  * of parameters.
  48.  */
  49.  
  50.  
  51. #define FIT_MAIN
  52.  
  53. #include <signal.h>
  54.  
  55. #include "plot.h"
  56. #include "matrix.h"
  57. #include "fit.h"
  58. #include "setshow.h"        /* for load_range */
  59. #include "alloc.h"
  60.  
  61. #if defined(MSDOS) || defined(DOS386)    /* non-blocking IO stuff */
  62. # include <io.h>
  63. # ifndef _Windows        /* WIN16 does define MSDOS .... */
  64. #  include <conio.h>
  65. # endif
  66. # include <dos.h>
  67. #else /* !(MSDOS || DOS386) */
  68. # ifndef VMS
  69. #  include <fcntl.h>
  70. # endif                /* !VMS */
  71. #endif /* !(MSDOS || DOS386) */
  72.  
  73. #if defined(ATARI) || defined(MTOS)
  74. # define getchx() Crawcin()
  75. static int kbhit(void);
  76. #endif
  77.  
  78. #define STANDARD    stderr    /* compatible with gnuplot philosophy */
  79.  
  80. #define BACKUP_SUFFIX ".old"
  81.  
  82.  
  83. /* access external global variables  (ought to make a globals.h someday) */
  84.  
  85. extern struct udft_entry *dummy_func;
  86. extern char dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
  87. extern char c_dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
  88. extern int c_token;
  89. extern int df_datum, df_line_number;
  90.  
  91. /* following 2 external arrays are needed to use time data */
  92. extern int datatype[];
  93. extern int df_timecol[];
  94.  
  95. enum marq_res {
  96.     OK, ERROR, BETTER, WORSE
  97. };
  98. typedef enum marq_res marq_res_t;
  99.  
  100. #ifdef INFINITY
  101. # undef INFINITY
  102. #endif
  103.  
  104. #define INFINITY    1e30
  105. #define NEARLY_ZERO 1e-30
  106.  
  107. /* create new variables with this value (was NEARLY_ZERO) */
  108. #define INITIAL_VALUE 1.0
  109.  
  110. /* Relative change for derivatives */
  111. #define DELTA        0.001
  112.  
  113. #define MAX_DATA    2048
  114. #define MAX_PARAMS  32
  115. #define MAX_LAMBDA  1e20
  116. #define MIN_LAMBDA  1e-20
  117. #define LAMBDA_UP_FACTOR 10
  118. #define LAMBDA_DOWN_FACTOR 10
  119.  
  120. #if defined(MSDOS) || defined(OS2) || defined(DOS386)
  121. # define PLUSMINUS   "\xF1"    /* plusminus sign */
  122. #else
  123. # define PLUSMINUS   "+/-"
  124. #endif
  125.  
  126. #define LASTFITCMDLENGTH 511
  127.  
  128. /* HBB 971023: new, allow for dynamic adjustment of these: */
  129. static int max_data;
  130. static int max_params;
  131.  
  132. static double epsilon = 1e-5;    /* convergence limit */
  133. static int maxiter = 0;        /* HBB 970304: maxiter patch */
  134.  
  135. static char fit_script[128];
  136. static char logfile[128] = "fit.log";
  137. static char *FIXED = "# FIXED";
  138. static char *GNUFITLOG = "FIT_LOG";
  139. static char *FITLIMIT = "FIT_LIMIT";
  140. static char *FITSTARTLAMBDA = "FIT_START_LAMBDA";
  141. static char *FITLAMBDAFACTOR = "FIT_LAMBDA_FACTOR";
  142. static char *FITMAXITER = "FIT_MAXITER";    /* HBB 970304: maxiter patch */
  143. static char *FITSCRIPT = "FIT_SCRIPT";
  144. static char *DEFAULT_CMD = "replot";    /* if no fitscript spec. */
  145. static char last_fit_command[LASTFITCMDLENGTH+1] = "";
  146.  
  147.  
  148. static FILE *log_f = NULL;
  149.  
  150. static int num_data, num_params;
  151. static int columns;
  152. static double *fit_x = 0, *fit_y = 0, *fit_z = 0, *err_data = 0, *a = 0;
  153. static TBOOLEAN ctrlc_flag = FALSE;
  154. static TBOOLEAN user_stop = FALSE;
  155.  
  156. static struct udft_entry func;
  157.  
  158. typedef char fixstr[MAX_ID_LEN+1];
  159.  
  160. static fixstr *par_name;
  161.  
  162. static double startup_lambda = 0, lambda_down_factor = LAMBDA_DOWN_FACTOR,
  163.  lambda_up_factor = LAMBDA_UP_FACTOR;
  164.  
  165. /*****************************************************************
  166.              internal Prototypes
  167. *****************************************************************/
  168.  
  169. static RETSIGTYPE ctrlc_handle __PROTO((int an_int));
  170. static void ctrlc_setup __PROTO((void));
  171. static void printmatrix __PROTO((double **C, int m, int n));
  172. static void print_matrix_and_vectors __PROTO((double **C, double *d, double *r, int m, int n));
  173. static marq_res_t marquardt __PROTO((double a[], double **alpha, double *chisq,
  174.                      double *lambda));
  175. static TBOOLEAN analyze __PROTO((double a[], double **alpha, double beta[],
  176.                  double *chisq));
  177. static void calculate __PROTO((double *zfunc, double **dzda, double a[]));
  178. static void call_gnuplot __PROTO((double *par, double *data));
  179. static TBOOLEAN fit_interrupt __PROTO((void));
  180. static TBOOLEAN regress __PROTO((double a[]));
  181. static void show_fit __PROTO((int i, double chisq, double last_chisq, double *a,
  182.                   double lambda, FILE * device));
  183. static TBOOLEAN is_empty __PROTO((char *s));
  184. static TBOOLEAN is_variable __PROTO((char *s));
  185. static double getdvar __PROTO((char *varname));
  186. static double createdvar __PROTO((char *varname, double value));
  187. static void splitpath __PROTO((char *s, char *p, char *f));
  188. static void backup_file __PROTO((char *, const char *));
  189.  
  190.  
  191. /*****************************************************************
  192.     Small function to write the last fit command into a file
  193.     Arg: Pointer to the file; if NULL, nothing is written,
  194.          but only the size of the string is returned.
  195. *****************************************************************/
  196.  
  197. size_t wri_to_fil_last_fit_cmd(fp)
  198. FILE *fp;
  199. {
  200.     if (fp == NULL)
  201.     return strlen(last_fit_command);
  202.     else
  203.     return (size_t) fputs(last_fit_command, fp);
  204. }
  205.  
  206.  
  207. /*****************************************************************
  208.     This is called when a SIGINT occurs during fit
  209. *****************************************************************/
  210.  
  211. static RETSIGTYPE ctrlc_handle(an_int)
  212. int an_int;
  213. {
  214. #ifdef OS2
  215.     (void) signal(an_int, SIG_ACK);
  216. #else
  217.     /* reinstall signal handler (necessary on SysV) */
  218.     (void) signal(SIGINT, (sigfunc) ctrlc_handle);
  219. #endif
  220.     ctrlc_flag = TRUE;
  221. }
  222.  
  223.  
  224. /*****************************************************************
  225.     setup the ctrl_c signal handler
  226. *****************************************************************/
  227. static void ctrlc_setup()
  228. {
  229. /*
  230.  *  MSDOS defines signal(SIGINT) but doesn't handle it through
  231.  *  real interrupts. So there remain cases in which a ctrl-c may
  232.  *  be uncaught by signal. We must use kbhit() instead that really
  233.  *  serves the keyboard interrupt (or write an own interrupt func
  234.  *  which also generates #ifdefs)
  235.  *
  236.  *  I hope that other OSes do it better, if not... add #ifdefs :-(
  237.  */
  238. #if (defined(__EMX__) || !defined(MSDOS) && !defined(DOS386)) && !defined(ATARI) && !defined(MTOS)
  239.     (void) signal(SIGINT, (sigfunc) ctrlc_handle);
  240. #endif
  241. }
  242.  
  243.  
  244. /*****************************************************************
  245.     getch that handles also function keys etc.
  246. *****************************************************************/
  247. #if defined(MSDOS) || defined(DOS386)
  248.  
  249. /* HBB 980317: added a prototype... */
  250. int getchx __PROTO((void));
  251.  
  252. int getchx()
  253. {
  254.     register int c = getch();
  255.     if (!c || c == 0xE0) {
  256.     c <<= 8;
  257.     c |= getch();
  258.     }
  259.     return c;
  260. }
  261. #endif
  262.  
  263. /*****************************************************************
  264.     in case of fatal errors
  265. *****************************************************************/
  266. void error_ex()
  267. {
  268.     char *sp;
  269.  
  270.     strncpy(fitbuf, "         ", 9);        /* start after GNUPLOT> */
  271.     sp = strchr(fitbuf, NUL);
  272.     while (*--sp == '\n')
  273.     ;
  274.     strcpy(sp + 1, "\n\n");    /* terminate with exactly 2 newlines */
  275.     fputs(fitbuf, STANDARD);
  276.     if (log_f) {
  277.     fprintf(log_f, "BREAK: %s", fitbuf);
  278.     (void) fclose(log_f);
  279.     log_f = NULL;
  280.     }
  281.     if (func.at) {
  282.     free(func.at);        /* release perm. action table */
  283.     func.at = (struct at_type *) NULL;
  284.     }
  285.     /* restore original SIGINT function */
  286. #if !defined(ATARI) && !defined(MTOS)
  287.     interrupt_setup();
  288. #endif
  289.  
  290.     bail_to_command_line();
  291. }
  292.  
  293.  
  294. /*****************************************************************
  295.     New utility routine: print a matrix (for debugging the alg.)
  296. *****************************************************************/
  297. static void printmatrix(C, m, n)
  298. double **C;
  299. int m, n;
  300. {
  301.     int i, j;
  302.  
  303.     for (i = 0; i < m; i++) {
  304.     for (j = 0; j < n - 1; j++)
  305.         Dblf2("%.8g |", C[i][j]);
  306.     Dblf2("%.8g\n", C[i][j]);
  307.     }
  308.     Dblf("\n");
  309. }
  310.  
  311. /**************************************************************************
  312.     Yet another debugging aid: print matrix, with diff. and residue vector
  313. **************************************************************************/
  314. static void print_matrix_and_vectors(C, d, r, m, n)
  315. double **C;
  316. double *d, *r;
  317. int m, n;
  318. {
  319.     int i, j;
  320.  
  321.     for (i = 0; i < m; i++) {
  322.     for (j = 0; j < n; j++)
  323.         Dblf2("%8g ", C[i][j]);
  324.     Dblf3("| %8g | %8g\n", d[i], r[i]);
  325.     }
  326.     Dblf("\n");
  327. }
  328.  
  329.  
  330. /*****************************************************************
  331.     Marquardt's nonlinear least squares fit
  332. *****************************************************************/
  333. static marq_res_t marquardt(a, C, chisq, lambda)
  334. double a[];
  335. double **C;
  336. double *chisq;
  337. double *lambda;
  338. {
  339.     int i, j;
  340.     static double *da = 0,    /* delta-step of the parameter */
  341.     *temp_a = 0,        /* temptative new params set   */
  342.     *d = 0, *tmp_d = 0, **tmp_C = 0, *residues = 0;
  343.     double tmp_chisq;
  344.  
  345.     /* Initialization when lambda == -1 */
  346.  
  347.     if (*lambda == -1) {    /* Get first chi-square check */
  348.     TBOOLEAN analyze_ret;
  349.  
  350.     temp_a = vec(num_params);
  351.     d = vec(num_data + num_params);
  352.     tmp_d = vec(num_data + num_params);
  353.     da = vec(num_params);
  354.     residues = vec(num_data + num_params);
  355.     tmp_C = matr(num_data + num_params, num_params);
  356.  
  357.     analyze_ret = analyze(a, C, d, chisq);
  358.  
  359.     /* Calculate a useful startup value for lambda, as given by Schwarz */
  360.     /* FIXME: this is doesn't turn out to be much better, really... */
  361.     if (startup_lambda != 0)
  362.         *lambda = startup_lambda;
  363.     else {
  364.         *lambda = 0;
  365.         for (i = 0; i < num_data; i++)
  366.         for (j = 0; j < num_params; j++)
  367.             *lambda += C[i][j] * C[i][j];
  368.         *lambda = sqrt(*lambda / num_data / num_params);
  369.     }
  370.  
  371.     /* Fill in the lower square part of C (the diagonal is filled in on
  372.        each iteration, see below) */
  373.     for (i = 0; i < num_params; i++)
  374.         for (j = 0; j < i; j++)
  375.         C[num_data + i][j] = 0, C[num_data + j][i] = 0;
  376.     /*printmatrix(C, num_data+num_params, num_params); */
  377.     return analyze_ret ? OK : ERROR;
  378.     }
  379.     /* once converged, free dynamic allocated vars */
  380.  
  381.     if (*lambda == -2) {
  382.     free(d);
  383.     free(tmp_d);
  384.     free(da);
  385.     free(temp_a);
  386.     free(residues);
  387.     free_matr(tmp_C);
  388.     return OK;
  389.     }
  390.     /* Givens calculates in-place, so make working copies of C and d */
  391.  
  392.     for (j = 0; j < num_data + num_params; j++)
  393.     memcpy(tmp_C[j], C[j], num_params * sizeof(double));
  394.     memcpy(tmp_d, d, num_data * sizeof(double));
  395.  
  396.     /* fill in additional parts of tmp_C, tmp_d */
  397.  
  398.     for (i = 0; i < num_params; i++) {
  399.     /* fill in low diag. of tmp_C ... */
  400.     tmp_C[num_data + i][i] = *lambda;
  401.     /* ... and low part of tmp_d */
  402.     tmp_d[num_data + i] = 0;
  403.     }
  404.     /* printmatrix(tmp_C, num_data+num_params, num_params); */
  405.  
  406.     /* FIXME: residues[] isn't used at all. Why? Should it be used? */
  407.  
  408.     Givens(tmp_C, tmp_d, da, residues, num_params + num_data, num_params, 1);
  409.     /*print_matrix_and_vectors (tmp_C, tmp_d, residues,
  410.        num_params+num_data, num_params); */
  411.  
  412.     /* check if trial did ameliorate sum of squares */
  413.  
  414.     for (j = 0; j < num_params; j++)
  415.     temp_a[j] = a[j] + da[j];
  416.  
  417.     if (!analyze(temp_a, tmp_C, tmp_d, &tmp_chisq)) {
  418.     /* FIXME: will never be reached: always returns TRUE */
  419.     return ERROR;
  420.     }
  421.  
  422.     if (tmp_chisq < *chisq) {    /* Success, accept new solution */
  423.     if (*lambda > MIN_LAMBDA) {
  424.         (void) putc('/', stderr);
  425.         *lambda /= lambda_down_factor;
  426.     }
  427.     *chisq = tmp_chisq;
  428.     for (j = 0; j < num_data; j++) {
  429.         memcpy(C[j], tmp_C[j], num_params * sizeof(double));
  430.         d[j] = tmp_d[j];
  431.     }
  432.     for (j = 0; j < num_params; j++)
  433.         a[j] = temp_a[j];
  434.     return BETTER;
  435.     } else {            /* failure, increase lambda and return */
  436.     (void) putc('*', stderr);
  437.     *lambda *= lambda_up_factor;
  438.     return WORSE;
  439.     }
  440. }
  441.  
  442.  
  443. /* FIXME: in the new code, this function doesn't really do enough to be
  444.  * useful. Maybe it ought to be deleted, i.e. integrated with
  445.  * calculate() ?
  446.  */
  447. /*****************************************************************
  448.     compute chi-square and numeric derivations
  449. *****************************************************************/
  450. static TBOOLEAN analyze(a, C, d, chisq)
  451. double a[];
  452. double **C;
  453. double d[];
  454. double *chisq;
  455. {
  456. /*
  457.  *  used by marquardt to evaluate the linearized fitting matrix C
  458.  *  and vector d, fills in only the top part of C and d
  459.  *  I don't use a temporary array zfunc[] any more. Just use
  460.  *  d[] instead.
  461.  */
  462.     int i, j;
  463.  
  464.     *chisq = 0;
  465.     calculate(d, C, a);
  466.  
  467.     for (i = 0; i < num_data; i++) {
  468.     /* note: order reversed, as used by Schwarz */
  469.     d[i] = (d[i] - fit_z[i]) / err_data[i];
  470.     *chisq += d[i] * d[i];
  471.     for (j = 0; j < num_params; j++)
  472.         C[i][j] /= err_data[i];
  473.     }
  474.     /* FIXME: why return a value that is always TRUE ? */
  475.     return TRUE;
  476. }
  477.  
  478.  
  479. /* To use the more exact, but slower two-side formula, activate the
  480.    following line: */
  481. /*#define TWO_SIDE_DIFFERENTIATION */
  482. /*****************************************************************
  483.     compute function values and partial derivatives of chi-square
  484. *****************************************************************/
  485. static void calculate(zfunc, dzda, a)
  486. double *zfunc;
  487. double **dzda;
  488. double a[];
  489. {
  490.     int k, p;
  491.     double tmp_a;
  492.     double *tmp_high, *tmp_pars;
  493. #ifdef TWO_SIDE_DIFFERENTIATION
  494.     double *tmp_low;
  495. #endif
  496.  
  497.     tmp_high = vec(num_data);    /* numeric derivations */
  498. #ifdef TWO_SIDE_DIFFERENTIATION
  499.     tmp_low = vec(num_data);
  500. #endif
  501.     tmp_pars = vec(num_params);
  502.  
  503.     /* first function values */
  504.  
  505.     call_gnuplot(a, zfunc);
  506.  
  507.     /* then derivatives */
  508.  
  509.     for (p = 0; p < num_params; p++)
  510.     tmp_pars[p] = a[p];
  511.     for (p = 0; p < num_params; p++) {
  512.     tmp_a = fabs(a[p]) < NEARLY_ZERO ? NEARLY_ZERO : a[p];
  513.     tmp_pars[p] = tmp_a * (1 + DELTA);
  514.     call_gnuplot(tmp_pars, tmp_high);
  515. #ifdef TWO_SIDE_DIFFERENTIATION
  516.     tmp_pars[p] = tmp_a * (1 - DELTA);
  517.     call_gnuplot(tmp_pars, tmp_low);
  518. #endif
  519.     for (k = 0; k < num_data; k++)
  520. #ifdef TWO_SIDE_DIFFERENTIATION
  521.         dzda[k][p] = (tmp_high[k] - tmp_low[k]) / (2 * tmp_a * DELTA);
  522. #else
  523.         dzda[k][p] = (tmp_high[k] - zfunc[k]) / (tmp_a * DELTA);
  524. #endif
  525.     tmp_pars[p] = a[p];
  526.     }
  527.  
  528. #ifdef TWO_SIDE_DIFFERENTIATION
  529.     free(tmp_low);
  530. #endif
  531.     free(tmp_high);
  532.     free(tmp_pars);
  533. }
  534.  
  535.  
  536. /*****************************************************************
  537.     call internal gnuplot functions
  538. *****************************************************************/
  539. static void call_gnuplot(par, data)
  540. double *par;
  541. double *data;
  542. {
  543.     int i;
  544.     struct value v;
  545.  
  546.     /* set parameters first */
  547.     for (i = 0; i < num_params; i++) {
  548.     (void) Gcomplex(&v, par[i], 0.0);
  549.     setvar(par_name[i], v);
  550.     }
  551.  
  552.     for (i = 0; i < num_data; i++) {
  553.     /* calculate fit-function value */
  554.     (void) Gcomplex(&func.dummy_values[0], fit_x[i], 0.0);
  555.     (void) Gcomplex(&func.dummy_values[1], fit_y[i], 0.0);
  556.     evaluate_at(func.at, &v);
  557.     if (undefined)
  558.         Eex("Undefined value during function evaluation");
  559.     data[i] = real(&v);
  560.     }
  561. }
  562.  
  563.  
  564. /*****************************************************************
  565.     handle user interrupts during fit
  566. *****************************************************************/
  567. static TBOOLEAN fit_interrupt()
  568. {
  569.     while (TRUE) {
  570.     fputs("\n\n(S)top fit, (C)ontinue, (E)xecute FIT_SCRIPT:  ", STANDARD);
  571.     switch (getc(stdin)) {
  572.  
  573.     case EOF:
  574.     case 's':
  575.     case 'S':
  576.         fputs("Stop.", STANDARD);
  577.         user_stop = TRUE;
  578.         return FALSE;
  579.  
  580.     case 'c':
  581.     case 'C':
  582.         fputs("Continue.", STANDARD);
  583.         return TRUE;
  584.  
  585.     case 'e':
  586.     case 'E':{
  587.         int i;
  588.         struct value v;
  589.         char *tmp;
  590.  
  591.         tmp = (fit_script != 0 && *fit_script) ? fit_script : DEFAULT_CMD;
  592.         fprintf(STANDARD, "executing: %s", tmp);
  593.         /* set parameters visible to gnuplot */
  594.         for (i = 0; i < num_params; i++) {
  595.             (void) Gcomplex(&v, a[i], 0.0);
  596.             setvar(par_name[i], v);
  597.         }
  598.         sprintf(input_line, tmp);
  599.         (void) do_line();
  600.         }
  601.     }
  602.     }
  603.     return TRUE;
  604. }
  605.  
  606.  
  607. /*****************************************************************
  608.     frame routine for the marquardt-fit
  609. *****************************************************************/
  610. static TBOOLEAN regress(a)
  611. double a[];
  612. {
  613.     double **covar, *dpar, **C, chisq, last_chisq, lambda;
  614.     int iter, i, j;
  615.     marq_res_t res;
  616.  
  617.     chisq = last_chisq = INFINITY;
  618.     C = matr(num_data + num_params, num_params);
  619.     lambda = -1;        /* use sign as flag */
  620.     iter = 0;            /* iteration counter  */
  621.  
  622.     /* ctrlc now serves as Hotkey */
  623.     ctrlc_setup();
  624.  
  625.     /* Initialize internal variables and 1st chi-square check */
  626.     if ((res = marquardt(a, C, &chisq, &lambda)) == ERROR)
  627.     Eex("FIT: error occured during fit");
  628.     res = BETTER;
  629.  
  630.     show_fit(iter, chisq, chisq, a, lambda, STANDARD);
  631.     show_fit(iter, chisq, chisq, a, lambda, log_f);
  632.  
  633.     /* MAIN FIT LOOP: do the regression iteration */
  634.  
  635.     /* HBB 981118: initialize new variable 'user_break' */
  636.     user_stop = FALSE;
  637.  
  638.     do {
  639. /*
  640.  *  MSDOS defines signal(SIGINT) but doesn't handle it through
  641.  *  real interrupts. So there remain cases in which a ctrl-c may
  642.  *  be uncaught by signal. We must use kbhit() instead that really
  643.  *  serves the keyboard interrupt (or write an own interrupt func
  644.  *  which also generates #ifdefs)
  645.  *
  646.  *  I hope that other OSes do it better, if not... add #ifdefs :-(
  647.  *  EMX does not have kbhit.
  648.  * 
  649.  *  HBB: I think this can be enabled for DJGPP V2. SIGINT is actually
  650.  *  handled there, AFAIK.
  651.  */
  652. #if ((defined(MSDOS) || defined(DOS386)) && !defined(__EMX__)) || defined(ATARI) || defined(MTOS)
  653.     if (kbhit()) {
  654.         do {
  655.         getchx();
  656.         } while (kbhit());
  657.         ctrlc_flag = TRUE;
  658.     }
  659. #endif
  660.  
  661.     if (ctrlc_flag) {
  662.         show_fit(iter, chisq, last_chisq, a, lambda, STANDARD);
  663.         ctrlc_flag = FALSE;
  664.         if (!fit_interrupt())    /* handle keys */
  665.         break;
  666.     }
  667.     if (res == BETTER) {
  668.         iter++;
  669.         last_chisq = chisq;
  670.     }
  671.     if ((res = marquardt(a, C, &chisq, &lambda)) == BETTER)
  672.         show_fit(iter, chisq, last_chisq, a, lambda, STANDARD);
  673.     } while ((res != ERROR)
  674.          && (lambda < MAX_LAMBDA)
  675.          && ((maxiter == 0) || (iter <= maxiter))
  676.          && (res == WORSE
  677.          || ((chisq > NEARLY_ZERO)
  678.              ? ((last_chisq - chisq) / chisq)
  679.              : (last_chisq - chisq)) > epsilon
  680.          )
  681.     );
  682.  
  683.     /* fit done */
  684.  
  685. #if !defined(ATARI) && !defined(MTOS)
  686.     /* restore original SIGINT function */
  687.     interrupt_setup();
  688. #endif
  689.  
  690.     /* HBB 970304: the maxiter patch: */
  691.     if ((maxiter > 0) && (iter > maxiter)) {
  692.     Dblf2("\nMaximum iteration count (%d) reached. Fit stopped.\n", maxiter);
  693.     } else if (user_stop) {
  694.     Dblf2("\nThe fit was stopped by the user after %d iterations.\n", iter);
  695.     } else {
  696.     Dblf2("\nAfter %d iterations the fit converged.\n", iter);
  697.     }
  698.  
  699.     Dblf2("final sum of squares of residuals : %g\n", chisq);
  700.     if (chisq > NEARLY_ZERO) {
  701.     Dblf2("rel. change during last iteration : %g\n\n", (chisq - last_chisq) / chisq);
  702.     } else {
  703.     Dblf2("abs. change during last iteration : %g\n\n", (chisq - last_chisq));
  704.     }
  705.  
  706.     if (res == ERROR)
  707.     Eex("FIT: error occured during fit");
  708.  
  709.     /* compute errors in the parameters */
  710.  
  711.     if (num_data == num_params) {
  712.     int i;
  713.  
  714.     Dblf("\nExactly as many data points as there are parameters.\n");
  715.     Dblf("In this degenerate case, all errors are zero by definition.\n\n");
  716.     Dblf("Final set of parameters \n");
  717.     Dblf("======================= \n\n");
  718.     for (i = 0; i < num_params; i++)
  719.         Dblf3("%-15.15s = %-15g\n", par_name[i], a[i]);
  720.     } else if (chisq < NEARLY_ZERO) {
  721.     int i;
  722.  
  723.     Dblf("\nHmmmm.... Sum of squared residuals is zero. Can't compute errors.\n\n");
  724.     Dblf("Final set of parameters \n");
  725.     Dblf("======================= \n\n");
  726.     for (i = 0; i < num_params; i++)
  727.         Dblf3("%-15.15s = %-15g\n", par_name[i], a[i]);
  728.     } else {
  729.     Dblf2("degrees of freedom (ndf) : %d\n",  num_data - num_params);
  730.     Dblf2("rms of residuals      (stdfit) = sqrt(WSSR/ndf)      : %g\n", sqrt(chisq / (num_data - num_params)));
  731.      Dblf2("variance of residuals (reduced chisquare) = WSSR/ndf : %g\n\n", chisq / (num_data - num_params));
  732.  
  733.     /* get covariance-, Korrelations- and Kurvature-Matrix */
  734.     /* and errors in the parameters                     */
  735.  
  736.     /* compute covar[][] directly from C */
  737.     Givens(C, 0, 0, 0, num_data, num_params, 0);
  738.     /*printmatrix(C, num_params, num_params); */
  739.  
  740.     /* Use lower square of C for covar */
  741.     covar = C + num_data;
  742.     Invert_RtR(C, covar, num_params);
  743.     /*printmatrix(covar, num_params, num_params); */
  744.  
  745.     /* calculate unscaled parameter errors in dpar[]: */
  746.     dpar = vec(num_params);
  747.     for (i = 0; i < num_params; i++) {
  748.         /* FIXME: can this still happen ? */
  749.         if (covar[i][i] <= 0.0)    /* HBB: prevent floating point exception later on */
  750.         Eex("Calculation error: non-positive diagonal element in covar. matrix");
  751.         dpar[i] = sqrt(covar[i][i]);
  752.     }
  753.  
  754.     /* transform covariances into correlations */
  755.     for (i = 0; i < num_params; i++) {
  756.         /* only lower triangle needs to be handled */
  757.         for (j = 0; j <= i; j++)
  758.         covar[i][j] /= dpar[i] * dpar[j];
  759.     }
  760.  
  761.     /* scale parameter errors based on chisq */
  762.     chisq = sqrt(chisq / (num_data - num_params));
  763.     for (i = 0; i < num_params; i++)
  764.         dpar[i] *= chisq;
  765.  
  766.     Dblf("Final set of parameters            Asymptotic Standard Error\n");
  767.     Dblf("=======================            ==========================\n\n");
  768.  
  769.     for (i = 0; i < num_params; i++) {
  770.         double temp =
  771.         (fabs(a[i]) < NEARLY_ZERO) ? 0.0 : fabs(100.0 * dpar[i] / a[i]);
  772.         Dblf6("%-15.15s = %-15g  %-3.3s %-12.4g (%.4g%%)\n",
  773.           par_name[i], a[i], PLUSMINUS, dpar[i], temp);
  774.     }
  775.  
  776.     Dblf("\n\ncorrelation matrix of the fit parameters:\n\n");
  777.     Dblf("               ");
  778.  
  779.     for (j = 0; j < num_params; j++)
  780.         Dblf2("%-6.6s ", par_name[j]);
  781.  
  782.     Dblf("\n");
  783.     for (i = 0; i < num_params; i++) {
  784.         Dblf2("%-15.15s", par_name[i]);
  785.         for (j = 0; j <= i; j++) {
  786.         /* Only print lower triangle of symmetric matrix */
  787.         Dblf2("%6.3f ", covar[i][j]);
  788.         }
  789.         Dblf("\n");
  790.     }
  791.  
  792.     free(dpar);
  793.     }
  794.  
  795.     /* HBB 990220: re-imported this snippet from older versions. Finally,
  796.      * some user noticed that it *is* necessary, after all. Not even
  797.      * Carsten Grammes himself remembered what it was for... :-(
  798.      * The thing is: the value of the last parameter is not reset to
  799.      * its original one after the derivatives have been calculated
  800.      */
  801.     /* restore last parameter's value (not done by calculate) */
  802.     {
  803.     struct value val;
  804.     Gcomplex (&val, a[num_params-1], 0.0);
  805.     setvar (par_name[num_params-1], val);
  806.      }
  807.  
  808.     /* call destructor for allocated vars */
  809.     lambda = -2;        /* flag value, meaning 'destruct!' */
  810.     (void) marquardt(a, C, &chisq, &lambda);
  811.  
  812.     free_matr(C);
  813.     return TRUE;
  814. }
  815.  
  816.  
  817. /*****************************************************************
  818.     display actual state of the fit
  819. *****************************************************************/
  820. static void show_fit(i, chisq, last_chisq, a, lambda, device)
  821. int i;
  822. double chisq;
  823. double last_chisq;
  824. double *a;
  825. double lambda;
  826. FILE *device;
  827. {
  828.     int k;
  829.  
  830.     fprintf(device, "\n\n\
  831. Iteration %d\n\
  832. WSSR        : %-15g   delta(WSSR)/WSSR   : %g\n\
  833. delta(WSSR) : %-15g   limit for stopping : %g\n\
  834. lambda      : %g\n\n%s parameter values\n\n",
  835.       i, chisq, chisq > NEARLY_ZERO ? (chisq - last_chisq) / chisq : 0.0,
  836.         chisq - last_chisq, epsilon, lambda,
  837.         (i > 0 ? "resultant" : "initial set of free"));
  838.     for (k = 0; k < num_params; k++)
  839.     fprintf(device, "%-15.15s = %g\n", par_name[k], a[k]);
  840. }
  841.  
  842.  
  843.  
  844. /*****************************************************************
  845.     is_empty: check for valid string entries
  846. *****************************************************************/
  847. static TBOOLEAN is_empty(s)
  848. char *s;
  849. {
  850.     while (*s == ' ' || *s == '\t' || *s == '\n')
  851.     s++;
  852.     return (TBOOLEAN) (*s == '#' || *s == '\0');
  853. }
  854.  
  855.  
  856. /*****************************************************************
  857.     get next word of a multi-word string, advance pointer
  858. *****************************************************************/
  859. char *get_next_word(s, subst)
  860. char **s;
  861. char *subst;
  862. {
  863.     char *tmp = *s;
  864.  
  865.     while (*tmp == ' ' || *tmp == '\t' || *tmp == '=')
  866.     tmp++;
  867.     if (*tmp == '\n' || *tmp == '\0')    /* not found */
  868.     return NULL;
  869.     if ((*s = strpbrk(tmp, " =\t\n")) == NULL)
  870.     *s = tmp + strlen(tmp);
  871.     *subst = **s;
  872.     *(*s)++ = '\0';
  873.     return tmp;
  874. }
  875.  
  876.  
  877. /*****************************************************************
  878.     check for variable identifiers
  879. *****************************************************************/
  880. static TBOOLEAN is_variable(s)
  881. char *s;
  882. {
  883.     while (*s != '\0') {
  884.     if (!isalnum((int) *s) && *s != '_')
  885.         return FALSE;
  886.     s++;
  887.     }
  888.     return TRUE;
  889. }
  890.  
  891. /*****************************************************************
  892.     first time settings
  893. *****************************************************************/
  894. void init_fit()
  895. {
  896.     func.at = (struct at_type *) NULL;    /* need to parse 1 time */
  897. }
  898.  
  899.  
  900. /*****************************************************************
  901.         Set a GNUPLOT user-defined variable
  902. ******************************************************************/
  903.  
  904. void setvar(varname, data)
  905. char *varname;
  906. struct value data;
  907. {
  908.     register struct udvt_entry *udv_ptr = first_udv, *last = first_udv;
  909.  
  910.     /* check if it's already in the table... */
  911.  
  912.     while (udv_ptr) {
  913.     last = udv_ptr;
  914.     if (!strcmp(varname, udv_ptr->udv_name))
  915.         break;
  916.     udv_ptr = udv_ptr->next_udv;
  917.     }
  918.  
  919.     if (!udv_ptr) {        /* generate new entry */
  920.     last->next_udv = udv_ptr = (struct udvt_entry *)
  921.         gp_alloc(sizeof(struct udvt_entry), "fit setvar");
  922.     udv_ptr->next_udv = NULL;
  923.     }
  924.     safe_strncpy(udv_ptr->udv_name, varname, sizeof(udv_ptr->udv_name));
  925.     udv_ptr->udv_value = data;
  926.     udv_ptr->udv_undef = FALSE;
  927. }
  928.  
  929.  
  930.  
  931. /*****************************************************************
  932.     Read INTGR Variable value, return 0 if undefined or wrong type
  933. *****************************************************************/
  934. int getivar(varname)
  935. char *varname;
  936. {
  937.     register struct udvt_entry *udv_ptr = first_udv;
  938.  
  939.     while (udv_ptr) {
  940.     if (!strcmp(varname, udv_ptr->udv_name))
  941.         return udv_ptr->udv_value.type == INTGR
  942.         ? udv_ptr->udv_value.v.int_val    /* valid */
  943.         : 0;        /* wrong type */
  944.     udv_ptr = udv_ptr->next_udv;
  945.     }
  946.     return 0;            /* not in table */
  947. }
  948.  
  949.  
  950. /*****************************************************************
  951.     Read DOUBLE Variable value, return 0 if undefined or wrong type
  952.    I dont think it's a problem that it's an integer - div
  953. *****************************************************************/
  954. static double getdvar(varname)
  955. char *varname;
  956. {
  957.     register struct udvt_entry *udv_ptr = first_udv;
  958.  
  959.     for (; udv_ptr; udv_ptr = udv_ptr->next_udv)
  960.     if (strcmp(varname, udv_ptr->udv_name) == 0)
  961.         return real(&(udv_ptr->udv_value));
  962.  
  963.     /* get here => not found */
  964.     return 0;
  965. }
  966.  
  967. /*****************************************************************
  968.    like getdvar, but
  969.    - convert it from integer to real if necessary
  970.    - create it with value INITIAL_VALUE if not found or undefined
  971. *****************************************************************/
  972. static double createdvar(varname, value)
  973. char *varname;
  974. double value;
  975. {
  976.     register struct udvt_entry *udv_ptr = first_udv;
  977.  
  978.     for (; udv_ptr; udv_ptr = udv_ptr->next_udv)
  979.     if (strcmp(varname, udv_ptr->udv_name) == 0) {
  980.         if (udv_ptr->udv_undef) {
  981.         udv_ptr->udv_undef = 0;
  982.         (void) Gcomplex(&udv_ptr->udv_value, value, 0.0);
  983.         } else if (udv_ptr->udv_value.type == INTGR) {
  984.         (void) Gcomplex(&udv_ptr->udv_value, (double) udv_ptr->udv_value.v.int_val, 0.0);
  985.         }
  986.         return real(&(udv_ptr->udv_value));
  987.     }
  988.     /* get here => not found */
  989.  
  990.     {
  991.     struct value tempval;
  992.     (void) Gcomplex(&tempval, value, 0.0);
  993.     setvar(varname, tempval);
  994.     }
  995.  
  996.     return value;
  997. }
  998.  
  999.  
  1000. /*****************************************************************
  1001.     Split Identifier into path and filename
  1002. *****************************************************************/
  1003. static void splitpath(s, p, f)
  1004. char *s;
  1005. char *p;
  1006. char *f;
  1007. {
  1008.     register char *tmp;
  1009.     tmp = s + strlen(s) - 1;
  1010.     while (*tmp != '\\' && *tmp != '/' && *tmp != ':' && tmp - s >= 0)
  1011.     tmp--;
  1012.     strcpy(f, tmp + 1);
  1013.     safe_strncpy(p, s, (size_t) (tmp - s + 1));
  1014.     p[tmp - s + 1] = NUL;
  1015. }
  1016.  
  1017.  
  1018. /*****************************************************************
  1019.     write the actual parameters to start parameter file
  1020. *****************************************************************/
  1021. void update(pfile, npfile)
  1022. char *pfile, *npfile;
  1023. {
  1024.     char fnam[256], path[256], sstr[256], pname[64], tail[127], *s = sstr,
  1025.     *tmp, c;
  1026.     char ifilename[256], *ofilename;
  1027.     FILE *of, *nf;
  1028.     double pval;
  1029.  
  1030.  
  1031.     /* update pfile npfile:
  1032.        if npfile is a valid file name,
  1033.        take pfile as input file and
  1034.        npfile as output file
  1035.      */
  1036.     if (VALID_FILENAME(npfile)) {
  1037.     safe_strncpy(ifilename, pfile, sizeof(ifilename));
  1038.     ofilename = npfile;
  1039.     } else {
  1040. #ifdef BACKUP_FILESYSTEM
  1041.     /* filesystem will keep original as previous version */
  1042.     safe_strncpy(ifilename, pfile, sizeof(ifilename));
  1043. #else
  1044.     backup_file(ifilename, pfile);    /* will Eex if it fails */
  1045. #endif
  1046.     ofilename = pfile;
  1047.     }
  1048.  
  1049.     /* split into path and filename */
  1050.     splitpath(ifilename, path, fnam);
  1051.     if (!(of = fopen(ifilename, "r")))
  1052.     Eex2("parameter file %s could not be read", ifilename);
  1053.  
  1054.     if (!(nf = fopen(ofilename, "w")))
  1055.     Eex2("new parameter file %s could not be created", ofilename);
  1056.  
  1057.     while (fgets(s = sstr, sizeof(sstr), of) != NULL) {
  1058.  
  1059.     if (is_empty(s)) {
  1060.         fputs(s, nf);    /* preserve comments */
  1061.         continue;
  1062.     }
  1063.     if ((tmp = strchr(s, '#')) != NULL) {
  1064.         safe_strncpy(tail, tmp, sizeof(tail));
  1065.         *tmp = NUL;
  1066.     } else
  1067.         strcpy(tail, "\n");
  1068.  
  1069.     tmp = get_next_word(&s, &c);
  1070.     if (!is_variable(tmp) || strlen(tmp) > MAX_ID_LEN) {
  1071.         (void) fclose(nf);
  1072.         (void) fclose(of);
  1073.         Eex2("syntax error in parameter file %s", fnam);
  1074.     }
  1075.     safe_strncpy(pname, tmp, sizeof(pname));
  1076.     /* next must be '=' */
  1077.     if (c != '=') {
  1078.         tmp = strchr(s, '=');
  1079.         if (tmp == NULL) {
  1080.         (void) fclose(nf);
  1081.         (void) fclose(of);
  1082.         Eex2("syntax error in parameter file %s", fnam);
  1083.         }
  1084.         s = tmp + 1;
  1085.     }
  1086.     tmp = get_next_word(&s, &c);
  1087.     if (!sscanf(tmp, "%lf", &pval)) {
  1088.         (void) fclose(nf);
  1089.         (void) fclose(of);
  1090.         Eex2("syntax error in parameter file %s", fnam);
  1091.     }
  1092.     if ((tmp = get_next_word(&s, &c)) != NULL) {
  1093.         (void) fclose(nf);
  1094.         (void) fclose(of);
  1095.         Eex2("syntax error in parameter file %s", fnam);
  1096.     }
  1097.     /* now modify */
  1098.  
  1099.     if ((pval = getdvar(pname)) == 0)
  1100.         pval = (double) getivar(pname);
  1101.  
  1102.     sprintf(sstr, "%g", pval);
  1103.     if (!strchr(sstr, '.') && !strchr(sstr, 'e'))
  1104.         strcat(sstr, ".0");    /* assure CMPLX-type */
  1105.  
  1106.     fprintf(nf, "%-15.15s = %-15.15s   %s", pname, sstr, tail);
  1107.     }
  1108.  
  1109.     if (fclose(nf) || fclose(of))
  1110.     Eex("I/O error during update");
  1111.  
  1112. }
  1113.  
  1114.  
  1115. /*****************************************************************
  1116.     Backup a file by renaming it to something useful. Return
  1117.     the new name in tofile
  1118. *****************************************************************/
  1119. static void backup_file(tofile, fromfile)
  1120. char *tofile;
  1121. const char *fromfile;
  1122. /* tofile must point to a char array[] or allocated data. See update() */
  1123. {
  1124. #if defined (WIN32) || defined(MSDOS) || defined(VMS)
  1125.     char *tmpn;
  1126. #endif
  1127.  
  1128.     /* win32 needs to have two attempts at the rename, since it may
  1129.      * be running on win32s with msdos 8.3 names
  1130.      */
  1131.  
  1132. /* first attempt, for all o/s other than MSDOS */
  1133.  
  1134. #ifndef MSDOS
  1135.  
  1136. 3Ge, fromfile);
  1137.  
  1138. #ifdef VMS
  1139.     /* replace all dots with _, since we will be adding the only
  1140.      * dot allowed in VMS names
  1141.      */
  1142.     while ((tmpn = strchr(tofile, '.')) != NULL)
  1143.     *tmpn = '_';
  1144. #endif /*VMS */
  1145.  
  1146.     strcat(tofile, BACKUP_SUFFIX);
  1147.  
  1148.     if (rename(fromfile, tofile) == 0)
  1149.     return;            /* hurrah */
  1150. #endif
  1151.  
  1152.  
  1153. #if defined (WIN32) || defined(MSDOS)
  1154.  
  1155.     /* first attempt for msdos. Second attempt for win32s */
  1156.  
  1157.     /* beware : strncpy is very dangerous since it does not guarantee
  1158.      * to terminate the string. Copy up to 8 characters. If exactly
  1159.      * chars were copied, the string is not terminated. If the
  1160.      * source string was shorter than 8 chars, no harm is done
  1161.      * (here) by writing to offset 8.
  1162.      */
  1163.     safe_strncpy(tofile, fromfile, 8);
  1164.     /* tofile[8] = NUL; */
  1165.  
  1166.     while ((tmpn = strchr(tofile, '.')) != NULL)
  1167.     *tmpn = '_';
  1168.  
  1169.     strcat(tofile, BACKUP_SUFFIX);
  1170.  
  1171.     if (rename(fromfile, tofile) == 0)
  1172.     return;            /* success */
  1173.  
  1174. #endif /* win32 || msdos */
  1175.  
  1176.     /* get here => rename failed. */
  1177.     Eex3("Could not rename file %s to %s", fromfile, tofile);
  1178. }
  1179.  
  1180.  
  1181. /*****************************************************************
  1182.     Interface to the classic gnuplot-software
  1183. *****************************************************************/
  1184.  
  1185. void do_fit()
  1186. {
  1187.     TBOOLEAN autorange_x = 3, autorange_y = 3;    /* yes */
  1188.     /* HBB 980401: new: z range specification */
  1189.     TBOOLEAN autorange_z = 3;
  1190.     double min_x, max_x;    /* range to fit */
  1191.     double min_y, max_y;    /* range to fit */
  1192.     /* HBB 980401: new: z range specification */
  1193.     double min_z, max_z;    /* range to fit */
  1194.     int dummy_x = -1, dummy_y = -1;    /* eg  fit [u=...] [v=...] */
  1195.     /* HBB 981210: memorize position of possible third [ : ] spec: */
  1196.     int zrange_token = -1;
  1197.  
  1198.     int i;
  1199.     double v[4];
  1200.     double tmpd;
  1201.     time_t timer;
  1202.     int token1, token2, token3;
  1203.     char *tmp;
  1204.  
  1205.     /* first look for a restricted x fit range... */
  1206.  
  1207.     if (equals(c_token, "[")) {
  1208.     c_token++;
  1209.     if (isletter(c_token)) {
  1210.         if (equals(c_token + 1, "=")) {
  1211.         dummy_x = c_token;
  1212.         c_token += 2;
  1213.         }
  1214.         /* else parse it as a xmin expression */
  1215.     }
  1216.     autorange_x = load_range(FIRST_X_AXIS, &min_x, &max_x, autorange_x);
  1217.     if (!equals(c_token, "]"))
  1218.         int_error("']' expected", c_token);
  1219.     c_token++;
  1220.     }
  1221.     /* ... and y */
  1222.  
  1223.     if (equals(c_token, "[")) {
  1224.     c_token++;
  1225.     if (isletter(c_token)) {
  1226.         if (equals(c_token + 1, "=")) {
  1227.         dummy_y = c_token;
  1228.         c_token += 2;
  1229.         }
  1230.         /* else parse it as a ymin expression */
  1231.     }
  1232.     autorange_y = load_range(FIRST_Y_AXIS, &min_y, &max_y, autorange_y);
  1233.     if (!equals(c_token, "]"))
  1234.         int_error("']' expected", c_token);
  1235.     c_token++;
  1236.     }
  1237.     /* HBB 980401: new: allow restricting the z range as well */
  1238.     /* ... and z */
  1239.  
  1240.     if (equals(c_token, "[")) {
  1241.     zrange_token = c_token++;
  1242.     if (isletter(c_token))
  1243.         /* do *not* allow the z range being given with a variable name */
  1244.         int_error("Can't re-name dependent variable", c_token);
  1245.     autorange_z = load_range(FIRST_Z_AXIS, &min_z, &max_z, autorange_z);
  1246.     if (!equals(c_token, "]"))
  1247.         int_error("']' expected", c_token);
  1248.     c_token++;
  1249. #if 0 /* HBB 981210: move this to a later point */
  1250.     } else {
  1251.     /* Just in case I muck up things below: make sure that the z
  1252.      * range is the same as the y range, if it didn't get specified
  1253.      */
  1254.     autorange_z = autorange_y;
  1255.     min_z = min_y;
  1256.     max_z = max_y;
  1257. #endif
  1258.     }
  1259.  
  1260.  
  1261.     /* now compile the function */
  1262.  
  1263.     token1 = c_token;
  1264.  
  1265.     if (func.at) {
  1266.     free(func.at);
  1267.     func.at = NULL;        /* in case perm_at() does int_error */
  1268.     }
  1269.     dummy_func = &func;
  1270.  
  1271.  
  1272.     /* set the dummy variable names */
  1273.  
  1274.     if (dummy_x >= 0)
  1275.     copy_str(c_dummy_var[0], dummy_x, MAX_ID_LEN);
  1276.     else
  1277.     strcpy(c_dummy_var[0], dummy_var[0]);
  1278.  
  1279.     if (dummy_y >= 0)
  1280.     copy_str(c_dummy_var[0], dummy_y, MAX_ID_LEN);
  1281.     else
  1282.     strcpy(c_dummy_var[1], dummy_var[1]);
  1283.  
  1284.     func.at = perm_at();
  1285.     dummy_func = NULL;
  1286.  
  1287.     token2 = c_token;
  1288.  
  1289.     /* use datafile module to parse the datafile and qualifiers */
  1290.  
  1291.     columns = df_open(4);    /* up to 4 using specs allowed */
  1292.     if (columns == 1)
  1293.     int_error("Need 2 to 4 using specs", c_token);
  1294.  
  1295.     /* The following patch was made by Remko Scharroo, 25-Mar-1999
  1296.      * We need to check if one of the columns is time data, like
  1297.      * in plot2d and plot3d */
  1298.  
  1299.     if (datatype[FIRST_X_AXIS] == TIME) {
  1300.     if (columns < 2)
  1301.         int_error("Need full using spec for x time data", c_token);
  1302.     df_timecol[0] = 1;
  1303.     }
  1304.     if (datatype[FIRST_Y_AXIS] == TIME) {
  1305.     if (columns < 1)
  1306.         int_error("Need using spec for y time data", c_token);
  1307.     df_timecol[1] = 1;
  1308.     }
  1309.     /* HBB 990326: added this check. Just in case some wants to fit
  1310.      * time/date data depending on two other variables ... */
  1311.     if (datatype[FIRST_Z_AXIS] == TIME) {
  1312.     if (columns < 4)
  1313.         int_error("Need full using spec for z time data", c_token);
  1314.     else
  1315.         df_timecol[2] = 1;
  1316.     }
  1317.     /* End of patch by Remko Scharroo */
  1318.  
  1319.     /* HBB 980401: if this is a single-variable fit, we shouldn't have
  1320.      * allowed a variable name specifier for 'y': */
  1321.     if ((dummy_y >= 0) && (columns < 4))
  1322.     int_error("Can't re-name 'y' in a one-variable fit", dummy_y);
  1323.  
  1324.     /* HBB 981210: two range specs mean different things, depending
  1325.      * on wether this is a 2D or 3D fit */
  1326.     if (columns < 4) {
  1327.     if (zrange_token != -1)
  1328.         int_error("Three range-specs not allowed in on-variable fit", zrange_token);
  1329.     else {
  1330.         /* 2D fit, 2 ranges: second range is for *z*, not y: */
  1331.         autorange_z = autorange_y;
  1332.         if (autorange_y & 1)
  1333.         min_z = min_y;
  1334.         if (autorange_y & 2)
  1335.         max_z = max_y;
  1336.     }
  1337.     }
  1338.  
  1339.     /* defer actually reading the data until we have parsed the rest
  1340.      * of the line */
  1341.  
  1342.     token3 = c_token;
  1343.  
  1344.     tmpd = getdvar(FITLIMIT);    /* get epsilon if given explicitly */
  1345.     if (tmpd < 1.0 && tmpd > 0.0)
  1346.     epsilon = tmpd;
  1347.  
  1348.     /* HBB 970304: maxiter patch */
  1349.     maxiter = getivar(FITMAXITER);
  1350.  
  1351.     /* get startup value for lambda, if given */
  1352.     tmpd = getdvar(FITSTARTLAMBDA);
  1353.  
  1354.     if (tmpd > 0.0) {
  1355.     startup_lambda = tmpd;
  1356.     printf("Lambda Start value set: %g\n", startup_lambda);
  1357.     }
  1358.  
  1359.     /* get lambda up/down factor, if given */
  1360.     tmpd = getdvar(FITLAMBDAFACTOR);
  1361.  
  1362.     if (tmpd > 0.0) {
  1363.     lambda_up_factor = lambda_down_factor = tmpd;
  1364.     printf("Lambda scaling factors reset:  %g\n", lambda_up_factor);
  1365.     }
  1366.     *fit_script = NUL;
  1367.     if ((tmp = getenv(FITSCRIPT)) != NULL)
  1368.     strcpy(fit_script, tmp);
  1369.  
  1370.     tmp = getenv(GNUFITLOG);    /* open logfile */
  1371.     if (tmp != NULL) {
  1372.     char *tmp2 = &tmp[strlen(tmp) - 1];
  1373.     if (*tmp2 == '/' || *tmp2 == '\\')
  1374.         sprintf(logfile, "%s%s", tmp, logfile);
  1375.     else
  1376.         strcpy(logfile, tmp);
  1377.     }
  1378.     if (!log_f && /* div */ !(log_f = fopen(logfile, "a")))
  1379.     Eex2("could not open log-file %s", logfile);
  1380.     fputs("\n\n*******************************************************************************\n", log_f);
  1381.     (void) time(&timer);
  1382.     fprintf(log_f, "%s\n\n", ctime(&timer));
  1383.     {
  1384.     char *line = NULL;
  1385.  
  1386.     m_capture(&line, token2, token3 - 1);
  1387.     fprintf(log_f, "FIT:    data read from %s\n", line);
  1388.     free(line);
  1389.     }
  1390.  
  1391.     max_data = MAX_DATA;
  1392.     fit_x = vec(max_data);    /* start with max. value */
  1393.     fit_y = vec(max_data);
  1394.     fit_z = vec(max_data);
  1395.  
  1396.     /* first read in experimental data */
  1397.  
  1398.     err_data = vec(max_data);
  1399.     num_data = 0;
  1400.  
  1401.     while ((i = df_readline(v, 4)) != EOF) {
  1402.     if (num_data >= max_data) {
  1403.         /* increase max_data by factor of 1.5 */
  1404.         max_data = (max_data * 3) / 2;
  1405.         if (0
  1406.         || !redim_vec(&fit_x, max_data)
  1407.         || !redim_vec(&fit_y, max_data)
  1408.         || !redim_vec(&fit_z, max_data)
  1409.         || !redim_vec(&err_data, max_data)
  1410.         ) {
  1411.         /* Some of the reallocations went bad: */
  1412.         df_close();
  1413.         Eex2("Out of memory in fit: too many datapoints (%d)?", max_data);
  1414.         } else {
  1415.         /* Just so we know that the routine is at work: */
  1416.         fprintf(STANDARD, "Max. number of data points scaled up to: %d\n",
  1417.             max_data);
  1418.         }
  1419.     }
  1420.     switch (i) {
  1421.     case DF_UNDEFINED:
  1422.     case DF_FIRST_BLANK:
  1423.     case DF_SECOND_BLANK:
  1424.         continue;
  1425.     case 0:
  1426.         Eex2("bad data on line %d of datafile", df_line_number);
  1427.         break;
  1428.     case 1:        /* only z provided */
  1429.         v[2] = v[0];
  1430.         v[0] = (double) df_datum;
  1431.         break;
  1432.     case 2:        /* x and z */
  1433.         v[2] = v[1];
  1434.         break;
  1435.  
  1436.         /* only if the explicitly asked for 4 columns do we
  1437.          * do a 3d fit. (We can get here if they didn't
  1438.          * specify a using spec, and the file has 4 columns
  1439.          */
  1440.     case 4:        /* x, y, z, error */
  1441.         if (columns == 4)
  1442.         break;
  1443.         /* else fall through */
  1444.     case 3:        /* x, z, error */
  1445.         v[3] = v[2];    /* error */
  1446.         v[2] = v[1];    /* z */
  1447.         break;
  1448.  
  1449.     }
  1450.  
  1451.     /* skip this point if it is out of range */
  1452.     if (!(autorange_x & 1) && (v[0] < min_x))
  1453.         continue;
  1454.     if (!(autorange_x & 2) && (v[0] > max_x))
  1455.         continue;
  1456.     /* HBB 980401: yrange is only relevant for 3d-fits! */
  1457.     if (columns == 4) {
  1458.         if (!(autorange_y & 1) && (v[1] < min_y))
  1459.         continue;
  1460.         if (!(autorange_y & 2) && (v[1] > max_y))
  1461.         continue;
  1462.     }
  1463.     /* HBB 980401: check *z* range for all fits */
  1464.     if (!(autorange_z & 1) && (v[2] < min_z))
  1465.         continue;
  1466.     if (!(autorange_z & 2) && (v[2] > max_z))
  1467.         continue;
  1468.  
  1469.     fit_x[num_data] = v[0];
  1470.     fit_y[num_data] = v[1];
  1471.     fit_z[num_data] = v[2];
  1472.  
  1473.     /* we only use error if _explicitly_ asked for by a
  1474.      * using spec
  1475.      */
  1476.     err_data[num_data++] = (columns > 2) ? v[3] : 1;
  1477.     }
  1478.     df_close();
  1479.  
  1480.     if (num_data <= 1)
  1481.     Eex("No data to fit ");
  1482.  
  1483.     /* now resize fields to actual length: */
  1484.     redim_vec(&fit_x, num_data);
  1485.     redim_vec(&fit_y, num_data);
  1486.     redim_vec(&fit_z, num_data);
  1487.     redim_vec(&err_data, num_data);
  1488.  
  1489.     fprintf(log_f, "        #datapoints = %d\n", num_data);
  1490.  
  1491.     if (columns < 3)
  1492.     fputs("        residuals are weighted equally (unit weight)\n\n", log_f);
  1493.  
  1494.     {
  1495.     char *line = NULL;
  1496.  
  1497.     m_capture(&line, token1, token2 - 1);
  1498.     fprintf(log_f, "function used for fitting: %s\n", line);
  1499.     free(line);
  1500.     }
  1501.  
  1502.     /* read in parameters */
  1503.  
  1504.     max_params = MAX_PARAMS;    /* HBB 971023: make this resizeable */
  1505.  
  1506.     if (!equals(c_token, "via"))
  1507.     int_error("Need via and either parameter list or file", c_token);
  1508.  
  1509.     a = vec(max_params);
  1510.     par_name = (fixstr *) gp_alloc((max_params + 1) * sizeof(fixstr), "fit param");
  1511.     num_params = 0;
  1512.  
  1513.     if (isstring(++c_token)) {    /* It's a parameter *file* */
  1514.     TBOOLEAN fixed;
  1515.     double tmp_par;
  1516.     char c, *s;
  1517.     char sstr[MAX_LINE_LEN+1];
  1518.     FILE *f;
  1519.  
  1520.     quote_str(sstr, c_token++, MAX_LINE_LEN);
  1521.  
  1522.     /* get parameters and values out of file and ignore fixed ones */
  1523.  
  1524.     fprintf(log_f, "fitted parameters and initial values from file: %s\n\n", sstr);
  1525.     if (!(f = fopen(sstr, "r")))
  1526.         Eex2("could not read parameter-file %s", sstr);
  1527.     while (TRUE) {
  1528.         if (!fgets(s = sstr, (int) sizeof(sstr), f))    /* EOF found */
  1529.         break;
  1530.         if ((tmp = strstr(s, FIXED)) != NULL) {    /* ignore fixed params */
  1531.         *tmp = NUL;
  1532.         fprintf(log_f, "FIXED:  %s\n", s);
  1533.         fixed = TRUE;
  1534.         } else
  1535.         fixed = FALSE;
  1536.         if ((tmp = strchr(s, '#')) != NULL)
  1537.         *tmp = NUL;
  1538.         if (is_empty(s))
  1539.         continue;
  1540.         tmp = get_next_word(&s, &c);
  1541.         if (!is_variable(tmp) || strlen(tmp) > MAX_ID_LEN) {
  1542.         (void) fclose(f);
  1543.         Eex("syntax error in parameter file");
  1544.         }
  1545.         safe_strncpy(par_name[num_params], tmp, sizeof(fixstr));
  1546.         /* next must be '=' */
  1547.         if (c != '=') {
  1548.         tmp = strchr(s, '=');
  1549.         if (tmp == NULL) {
  1550.             (void) fclose(f);
  1551.             Eex("syntax error in parameter file");
  1552.         }
  1553.         s = tmp + 1;
  1554.         }
  1555.         tmp = get_next_word(&s, &c);
  1556.         if (sscanf(tmp, "%lf", &tmp_par) != 1) {
  1557.         (void) fclose(f);
  1558.         Eex("syntax error in parameter file");
  1559.         }
  1560.         /* make fixed params visible to GNUPLOT */
  1561.         if (fixed) {
  1562.         struct value tempval;
  1563.         (void) Gcomplex(&tempval, tmp_par, 0.0);
  1564.         /* use parname as temp */
  1565.         setvar(par_name[num_params], tempval);
  1566.         } else {
  1567.         if (num_params >= max_params) {
  1568.             max_params = (max_params * 3) / 2;
  1569.             if (0
  1570.             || !redim_vec(&a, max_params)
  1571.             || !(par_name = gp_realloc(par_name, (max_params + 1) * sizeof(fixstr), "fit param resize"))
  1572.             ) {
  1573.             (void) fclose(f);
  1574.             Eex("Out of memory in fit: too many parameters?");
  1575.             }
  1576.         }
  1577.         a[num_params++] = tmp_par;
  1578.         }
  1579.  
  1580.         if ((tmp = get_next_word(&s, &c)) != NULL) {
  1581.         (void) fclose(f);
  1582.         Eex("syntax error in parameter file");
  1583.         }
  1584.     }
  1585.     (void) fclose(f);
  1586.  
  1587.     } else {
  1588.     /* not a string after via: it's a variable listing */
  1589.  
  1590.     fputs("fitted parameters initialized with current variable values\n\n", log_f);
  1591.     do {
  1592.         if (!isletter(c_token))
  1593.         Eex("no parameter specified");
  1594.         capture(par_name[num_params], c_token, c_token, (int) sizeof(par_name[0]));
  1595.         if (num_params >= max_params) {
  1596.         max_params = (max_params * 3) / 2;
  1597.         if (0
  1598.             || !redim_vec(&a, max_params)
  1599.             || !(par_name = gp_realloc(par_name, (max_params + 1) * sizeof(fixstr), "fit param resize"))
  1600.             ) {
  1601.             Eex("Out of memory in fit: too many parameters?");
  1602.         }
  1603.         }
  1604.         /* create variable if it doesn't exist */
  1605.         a[num_params] = createdvar(par_name[num_params], INITIAL_VALUE);
  1606.         ++num_params;
  1607.     } while (equals(++c_token, ",") && ++c_token);
  1608.     }
  1609.  
  1610.     redim_vec(&a, num_params);
  1611.     par_name = (fixstr *) gp_realloc(par_name, (num_params + 1) * sizeof(fixstr), "fit param");
  1612.  
  1613.     if (num_data < num_params)
  1614.     Eex("Number of data points smaller than number of parameters");
  1615.  
  1616.     /* avoid parameters being equal to zero */
  1617.     for (i = 0; i < num_params; i++)
  1618.     if (a[i] == 0)
  1619.         a[i] = NEARLY_ZERO;
  1620.  
  1621.     (void) regress(a);
  1622.  
  1623.     (void) fclose(log_f);
  1624.     log_f = NULL;
  1625.     free(fit_x);
  1626.     free(fit_y);
  1627.     free(fit_z);
  1628.     free(err_data);
  1629.     free(a);
  1630.     free(par_name);
  1631.     if (func.at) {
  1632.     free(func.at);        /* release perm. action table */
  1633.     func.at = (struct at_type *) NULL;
  1634.     }
  1635.     safe_strncpy(last_fit_command, input_line, sizeof(last_fit_command));
  1636. }
  1637.  
  1638. #if defined(ATARI) || defined(MTOS)
  1639. static int kbhit()
  1640. {
  1641.     fd_set rfds;
  1642.     struct timeval timeout;
  1643.  
  1644.     timeout.tv_sec = 0;
  1645.     timeout.tv_usec = 0;
  1646.     rfds = (fd_set) (1L << fileno(stdin));
  1647.     return ((select(0, &rfds, NULL, NULL, &timeout) > 0) ? 1 : 0);
  1648. }
  1649. #endif
  1650.